{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96f94593",
   "metadata": {},
   "outputs": [],
   "source": [
    "remote = \"https://raw.githubusercontent.com/nansencenter/DA-tutorials\"\n",
    "!wget -qO- {remote}/master/notebooks/resources/colab_bootstrap.sh | bash -s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cb7f04e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from resources import show_answer, interact\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy.random as rnd\n",
    "plt.ion();\n",
    "rnd.seed(3000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ff70496",
   "metadata": {},
   "source": [
    "# T2 - The Gaussian (Normal) distribution\n",
    "\n",
    "We begin by reviewing the most useful of probability distributions.\n",
    "But first, let's refresh some basic theory.\n",
    "$\n",
    "\\newcommand{\\Reals}{\\mathbb{R}}\n",
    "\\newcommand{\\Expect}[0]{\\mathbb{E}}\n",
    "\\newcommand{\\NormDist}{\\mathscr{N}}\n",
    "\\newcommand{\\mat}[1]{{\\mathbf{{#1}}}}\n",
    "\\newcommand{\\vect}[1]{{\\mathbf{#1}}}\n",
    "\\newcommand{\\trsign}{{\\mathsf{T}}}\n",
    "\\newcommand{\\tr}{^{\\trsign}}\n",
    "\\newcommand{\\z}[0]{\\vect{z}}\n",
    "\\newcommand{\\E}[0]{\\mat{E}}\n",
    "\\newcommand{\\I}[0]{\\mat{I}}\n",
    "\\newcommand{\\x}[0]{\\vect{x}}\n",
    "\\newcommand{\\X}[0]{\\mat{X}}\n",
    "$\n",
    "\n",
    "<a name=\"Probability-essentials\"></a>\n",
    "\n",
    "## Probability essentials\n",
    "\n",
    "As stated by James Bernoulli (1713) and elucidated by [Laplace (1812)](#References):\n",
    "\n",
    "> The Probability for an event is the ratio of the number of cases favorable to it, to the number of all\n",
    "> cases possible when nothing leads us to expect that any one of these cases should occur more than any other,\n",
    "> which renders them, for us, equally possible:\n",
    "\n",
    "$$ \\mathbb{P}(\\text{event}) = \\frac{\\text{number of} \\textit{ favorable } \\text{outcomes}}{\\text{number of} \\textit{ possible } \\text{outcomes}} $$\n",
    "\n",
    "The probability of *both* events $A$ and $B$ occurring is given by their intersection:\n",
    "$\\mathbb{P}(A \\cap B)$, while the probability of *either (or)* is obtained by their union $\\mathbb{P}(A \\cup B)$.\n",
    "The *conditional* probability of $A$ given $B$ restricts our attention (count)\n",
    "to cases where $B$ occurs: $\\mathbb{P}(A | B) = \\frac{\\mathbb{P}(A \\cap B)}{\\mathbb{P}(B)}$.\n",
    "\n",
    "A **random variable**, $X$, is a *numeric quantity* taking values as a function of some underlying random process.\n",
    "Rather than \"*did* event $A$ occur or no?\",\n",
    "random variables conveniently enable the question\n",
    "\"*what* was the value of $X$?\".\n",
    "Each value, $x$, constitutes an event that is disjoint from all others (functions never being one-to-many),\n",
    "and so they define a probability space of outcomes with associated probabilities,\n",
    "which can be tabulated into *distributions*.\n",
    "If $X$ is *discrete*, then $p_X(x) := \\mathbb{P}(X{=}x)$ is a list mapping outcomes to probabilities\n",
    "called the probability *mass* function (**pmf**).\n",
    "It sums to 1, and may be written $p(x)$ if contextually unambiguous.\n",
    "The cumulative distribution function (**cdf**) is defined as $F(x) := \\mathbb{P}(X \\le x)$.\n",
    "The 2D table of *joint* probabilities of $X$ and $Y$ is denoted $p(x, y) = \\mathbb{P}(X{=}x \\cap Y{=}y)$,\n",
    "while the conditionals are denoted $p(x|y) = \\frac{p(x,y)}{p(y)}$.\n",
    "\n",
    "- The *marginal* pmf, $p(x)$, can be recovered from the joint pmf, $p(x, y)$, by summing over all $y$.\n",
    "- *Independence* means $p(x, y) = p(x) \\, p(y)$ for all possible $x, y$. Equivalently $p(x|y) = p(x)$.\n",
    "\n",
    "We will mainly be concerned with *continuous* random variables,\n",
    "for which $\\mathbb{P}(X \\in I)$ may be non-zero for any interval, $I$.\n",
    "The distribution of $X$ is then characterised by its probability *density* function (**pdf**),\n",
    "defined as $p(x) = F'(x)$ or\n",
    "\n",
    "$$p(x) = \\lim_{h \\to 0} \\frac{\\mathbb{P}(X \\in [x,\\, x{+} h])}{h} \\,.$$\n",
    "\n",
    "The **sample average** of draws from a random variable $X$\n",
    "is denoted with an overhead bar:\n",
    "$$ \\bar{x} := \\frac{1}{N} \\sum_{n=1}^{N} x_n \\,. $$\n",
    "The *law of large numbers (LLN)* states that, as $N \\to \\infty$,\n",
    "the sample average converges to the **expected value** (sometimes called the **mean**):\n",
    "$$ \\Expect[X] ≔ \\int x \\, p(x) \\, d x \\,, $$\n",
    "where the (omitted) domain of integration is *all values of $x$*.\n",
    "\n",
    "## The univariate (a.k.a. 1-dimensional, scalar) Gaussian\n",
    "\n",
    "If $X$ is Gaussian (also known as \"Normal\"), we write\n",
    "$X \\sim \\NormDist(\\mu, \\sigma^2)$, or $p(x) = \\NormDist(x \\mid \\mu, \\sigma^2)$,\n",
    "where the parameters $\\mu$ and $\\sigma^2$ are called the mean and variance\n",
    "(for reasons that will become clear below).\n",
    "The Gaussian pdf, for $x \\in (-\\infty, +\\infty)$, is\n",
    "$$ \\large \\NormDist(x \\mid \\mu, \\sigma^2) = (2 \\pi \\sigma^2)^{-1/2} e^{-(x-\\mu)^2/2 \\sigma^2} \\, . \\tag{G1} $$\n",
    "\n",
    "Run the cell below to define a function to compute the pdf (G1) using the `scipy` library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3be2918",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pdf_G1(x, mu, sigma2):\n",
    "    \"Univariate Gaussian pdf\"\n",
    "    pdf_values = sp.stats.norm.pdf(x, loc=mu, scale=np.sqrt(sigma2))\n",
    "    return pdf_values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "069911f0",
   "metadata": {},
   "source": [
    "Computers typically represent functions *numerically* by their values at a set of grid points (nodes),\n",
    "an approach called ***discretization***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80586aed",
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = -20, 20\n",
    "N = 201                         # num of grid points\n",
    "grid1d = np.linspace(*bounds,N) # grid\n",
    "dx = grid1d[1] - grid1d[0]      # grid spacing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e40b3ce",
   "metadata": {},
   "source": [
    "Feel free to return here later and change the grid resolution to see how\n",
    "it affects the cells below (after re-running them).\n",
    "\n",
    "The following code plots the Gaussian pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0442e503",
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = []\n",
    "@interact(mu=bounds, sigma=(.1, 10, 1))\n",
    "def plot_pdf(mu=0, sigma=5):\n",
    "    plt.figure(figsize=(6, 2))\n",
    "    colors = plt.get_cmap('hsv')([(k-len(hist))%9/9 for k in range(9)])\n",
    "    plt.xlim(*bounds)\n",
    "    plt.ylim(0, .2)\n",
    "    hist.insert(0, pdf_G1(grid1d, mu, sigma**2))\n",
    "    for density_values, color in zip(hist, colors):\n",
    "        plt.plot(grid1d, density_values, c=color)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e0059c2",
   "metadata": {},
   "source": [
    "#### Exc – parameter influence\n",
    "\n",
    "Experiment with `mu` and `sigma` to answer these questions:\n",
    "\n",
    "- How does the pdf curve change when `mu` changes? (Several options may be correct or incorrect)\n",
    "<details style=\"border: 1px solid #aaaaaa; border-radius: 4px; padding: 0.5em 0.5em 0;\">\n",
    "<summary style=\"font-weight: normal; font-style: italic; margin: -0.5em -0.5em 0; padding: 0.5em;\">\n",
    "  Click to view options 🔍\n",
    "</summary>\n",
    "\n",
    "  1. It changes the curve into a uniform distribution.\n",
    "  1. It changes the width of the curve.\n",
    "  1. It shifts the peak of the curve to the left or right.\n",
    "  1. It changes the height of the curve.\n",
    "  1. It transforms the curve into a binomial distribution.\n",
    "  1. It makes the curve wider or narrower.\n",
    "  1. It modifies the skewness (asymmetry) of the curve.\n",
    "  1. It causes the curve to expand vertically while keeping the width the same.\n",
    "  1. It translates the curve horizontally.\n",
    "  1. It alters the kurtosis (peakedness) of the curve.\n",
    "  1. It rotates the curve around the origin.\n",
    "  1. It makes the curve a straight line.\n",
    "</details>\n",
    "\n",
    "- How does the pdf curve change when you increase `sigma`?  \n",
    "  Refer to the same options as the previous question.\n",
    "- In a few words, describe the shape of the Gaussian pdf curve.\n",
    "  Does this remind you of anything? *Hint: it should be clear as a bell!*\n",
    "\n",
    "**Exc – Implementation:** Change the implementation of `pdf_G1` so that it does not use `scipy`, but instead uses your own code (with `numpy` only). Re-run all of the above cells and check that you get the same plots as before.  \n",
    "*Hint: `**` is the exponentiation/power operator, but $e^x$ is more efficiently computed with `np.exp(x)`*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3e1106d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('pdf_G1')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94b6d541",
   "metadata": {},
   "source": [
    "**Exc – Derivatives:** Recall $p(x) = \\NormDist(x \\mid \\mu, \\sigma^2)$ from eqn. (G1).  \n",
    "Use pen, paper, and calculus to answer the following questions,\n",
    "which will help you remember some key properties of the distribution.\n",
    "\n",
    "- (i) Find $x$ such that $p(x) = 0$.\n",
    "- (ii) Where is the location of the **mode (maximum)** of the density?  \n",
    "  I.e. find $x$ such that $\\frac{d p}{d x}(x) = 0$.\n",
    "  *Hint: begin by writing $p(x)$ as $c e^{- J(x)}$ for some $J(x)$.*\n",
    "- (iii) Where is the **inflection point**? I.e. where $\\frac{d^2 p}{d x^2}(x) = 0$.\n",
    "- (iv) *Optional*: Some forms of *sensitivity analysis* (typically for non-Gaussian $p$) consist in estimating/approximating the Hessian, i.e. $\\frac{d^2 \\log p}{d x^2}$. Explain what this has to do with *uncertainty quantification*.\n",
    "\n",
    "<a name=\"Exc-(optional)-–-Change-of-variables\"></a>\n",
    "\n",
    "#### Exc (optional) – Change of variables\n",
    "\n",
    "Let $U = \\phi(X)$ for some monotonic function $\\phi$,\n",
    "and let $p_x$ and $p_u$ be their probability density functions (pdf).\n",
    "\n",
    "- (a): Show that $p_u(u) = p_x\\big(\\phi^{-1}(u)\\big) \\frac{1}{|\\phi'(u)|}$,\n",
    "- (b): Show that you don't need to derive the density of $u$ in order to compute its expectation, i.e. that\n",
    "  $$ \\Expect[U] = \\int  \\phi(x) \\, p_x(x) \\, d x ≕ \\Expect[\\phi(x)] \\,,$$\n",
    "  *PS: this result is [pretty intuitive](https://en.wikipedia.org/wiki/Law_of_the_unconscious_statistician),\n",
    "  and also holds for non-injective transformations, $\\phi$,\n",
    "  as well as functions of multiple random variables.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92777530",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('CVar in proba')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8311df4c",
   "metadata": {},
   "source": [
    "<a name=\"Exc-(optional)-–-Integrals\"></a>\n",
    "\n",
    "#### Exc (optional) – Integrals\n",
    "\n",
    "Recall $p(x) = \\NormDist(x \\mid \\mu, \\sigma^2)$ from eqn. (G1). Abbreviate it as $c = (2 \\pi \\sigma^2)^{-1/2}$.  \n",
    "Use pen, paper, and calculus to show that\n",
    "\n",
    "- (i) the first parameter, $\\mu$, indicates its **mean**, i.e. that $$\\mu = \\Expect[X] \\,.$$\n",
    "  *Hint: you can rely on the result of (iii)*\n",
    "- (ii) the second parameter, $\\sigma^2>0$, indicates its **variance**,\n",
    "  i.e. that $$\\sigma^2 = \\mathbb{Var}(X) \\mathrel{≔} \\Expect[(X-\\mu)^2] \\,.$$\n",
    "  *Hint: use $x^2 = x x$ to enable integration by parts.*\n",
    "- (iii) $c$ is indeed the right normalizing constant, i.e. that\n",
    "  $$E[1] = 1 \\,.$$\n",
    "  *Hint: Neither Bernoulli and Laplace managed this,\n",
    "  until [Gauss (1809)](#References) did by first deriving $(E[1])^2$.\n",
    "  Here is a nice [video demonstration by 3Blue1Brown](https://www.youtube.com/watch?v=cy8r7WSuT1I&t=3m52s).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5470edc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('Gauss integrals')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce40fa45",
   "metadata": {},
   "source": [
    "**Exc (optional) – Riemann sums**:\n",
    "Recall that integrals (for example for the mean and variance)\n",
    "compute an \"area under the curve\".\n",
    "On a discrete grid, integrals can be approximated using the [Trapezoidal rule](https://en.wikipedia.org/wiki/Riemann_sum#Trapezoidal_rule).\n",
    "\n",
    "- (a) Replace `np.trapezoid` below with your own implementation (using `sum()`).\n",
    "- (b) Use `np.trapezoid` to compute the probability that a scalar Gaussian $X$ lies within $1$ standard deviation of its mean.\n",
    "  *Hint: the numerical answer you should find is $\\mathbb{P}(X \\in [\\mu {-} \\sigma, \\mu {+} \\sigma]) \\approx 68\\%$.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05f83331",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_and_var(pdf_values, grid):\n",
    "    f, x = pdf_values, grid\n",
    "    mu = np.trapezoid(f*x, x)\n",
    "    s2 = np.trapezoid(f*(x-mu)**2, x)\n",
    "    return mu, s2\n",
    "\n",
    "mu, sigma = 0, 2 # example\n",
    "pdf_vals = pdf_G1(grid1d, mu=mu, sigma2=sigma**2)\n",
    "'Should equal mu and sigma2: %f, %f' % mean_and_var(pdf_vals, grid1d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cf84c4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('Riemann sums', 'a')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a46ec374",
   "metadata": {},
   "source": [
    "**Exc – The uniform pdf**:\n",
    "Below is the pdf of the [uniform/flat/box distribution](https://en.wikipedia.org/wiki/Uniform_distribution_(continuous))\n",
    "for a given mean and variance.\n",
    "\n",
    "- Use `mean_and_var()` to verify `pdf_U1` (as is).\n",
    "- Replace `_G1` with `_U1` in the code generating the above interactive plot.\n",
    "- Why are the walls (ever so slightly) inclined?\n",
    "- Write your own implementation below, and check that it reproduces the `scipy` version already in place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35c1e382",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pdf_U1(x, mu, sigma2):\n",
    "    a = mu - np.sqrt(3*sigma2)\n",
    "    b = mu + np.sqrt(3*sigma2)\n",
    "    pdf_values = sp.stats.uniform(loc=a, scale=(b-a)).pdf(x)\n",
    "    # Your own implementation:\n",
    "    # height = ...\n",
    "    # pdf_values = height * np.ones_like(x)\n",
    "    # pdf_values[x<a] = ...\n",
    "    # pdf_values[x>b] = ...\n",
    "    return pdf_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94e09bbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('pdf_U1')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b28223c0",
   "metadata": {},
   "source": [
    "## The multivariate (i.e. vector) Gaussian\n",
    "\n",
    "A *multivariate* random variable, i.e. a **vector**, is simply a collection of scalar variables (on the same probability space).\n",
    "Its distribution is the *joint* distribution of its components.\n",
    "The pdf of the multivariate Gaussian $\\X$ (for any dimension $\\ge 1$) is\n",
    "\n",
    "$$\\large \\NormDist(\\x \\mid \\mathbf{\\mu}, \\mathbf{\\Sigma}) =\n",
    "|2 \\pi \\mathbf{\\Sigma}|^{-1/2} \\, \\exp\\Big(-\\frac{1}{2}\\|\\x-\\mathbf{\\mu}\\|^2_\\mathbf{\\Sigma} \\Big) \\,, \\tag{GM} $$\n",
    "which is very similar to the univariate (scalar) case (G1),\n",
    "but with $|.|$ representing the matrix determinant,\n",
    "and $\\|.\\|_\\mathbf{W}$ representing the weighted 2-norm: $\\|\\x\\|^2_\\mathbf{W} = \\x^T \\mathbf{W}^{-1} \\x$.  \n",
    "\n",
    "<details style=\"border: 1px solid #aaaaaa; border-radius: 4px; padding: 0.5em 0.5em 0;\">\n",
    "<summary style=\"font-weight: normal; font-style: italic; margin: -0.5em -0.5em 0; padding: 0.5em;\">\n",
    "  $\\mathbf{W}$ must be symmetric-positive-definite (SPD) because ... (optional reading 🔍)\n",
    "</summary>\n",
    "\n",
    "- The norm (a quadratic form) is invariant to any asymmetry in the weight matrix.\n",
    "- The density (GM) would not be integrable (over $\\Reals^{d}$) if $\\x\\tr \\mathbf{\\Sigma} \\x > 0$.\n",
    "\n",
    "- - -\n",
    "</details>\n",
    "\n",
    "Moreover, [as above](#Exc-(optional)-–-Integrals), it can be shown that\n",
    "\n",
    "- $\\mathbf{\\mu} = \\Expect[\\X]$,\n",
    "- $\\mathbf{\\Sigma} = \\Expect[(\\X-\\mu)(\\X-\\mu)\\tr]  =: \\mathbb{Cov}(\\X)$.\n",
    "\n",
    "As such, $\\mathbf{\\Sigma}$ is called the **covariance matrix**,\n",
    "whose individual elements are individual covariances,\n",
    "$\\Sigma_{i,j} = \\Expect[(X_i-\\mu_i)(X_j-\\mu_j)] =: \\mathbb{Cov}(X_i, X_j)$,\n",
    "and – on the diagonal – variances: $\\Sigma_{i,i} = \\mathbb{Var}(X_i)$.\n",
    "\n",
    "The following implements the pdf (GM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f75d1bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.linalg as la\n",
    "\n",
    "def pdf_GM(x, mu, Sigma):\n",
    "    \"pdf – Gaussian, Multivariate: N(x | mu, Sigma) for each x.\"\n",
    "    c = np.sqrt(la.det(2*np.pi*Sigma))\n",
    "    return 1/c * np.exp(-0.5*weighted_norm22(x - mu, Sigma))\n",
    "\n",
    "def weighted_norm22(points, cov):\n",
    "    \"Computes the weighted norm of each vector (a row in `points`).\"\n",
    "    W = la.inv(cov) # NB: replace by la.solve() in real applications!\n",
    "    return np.sum( (points @ W) * points, axis=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74bb5cf0",
   "metadata": {},
   "source": [
    "The norm implementation is a bit tricky because it uses `@` (matrix multiplication), `*` (array, i.e. element-wise multiplication) and `axis=-1` (sum along the last dimension) to enable inputting the entire lattice/grid at once without shape manipulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "affc16da",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid2d = np.dstack(np.meshgrid(grid1d, grid1d))\n",
    "grid2d.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52569aaf",
   "metadata": {},
   "source": [
    "The following code plots the pdf as contour (level) curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dc581fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(corr=(-1, 1, .001), std_x=(1e-5, 10, 1), seed=(0, 9))\n",
    "def plot_pdf_G2(corr=0.7, std_x=1, seed=0):\n",
    "    mu = 0\n",
    "    var_x = std_x**2\n",
    "    var_y = 1\n",
    "    cv_xy = np.sqrt(var_x * var_y) * corr\n",
    "\n",
    "    # Assemble covariance matrix (C)\n",
    "    C = 25 * np.array([[var_x, cv_xy],\n",
    "                       [cv_xy, var_y]])\n",
    "\n",
    "    # Evaluate (compute)\n",
    "    density_values = pdf_GM(grid2d, mu=mu, Sigma=C)\n",
    "\n",
    "    # Plot\n",
    "    plt.figure(figsize=(4, 4))\n",
    "    plt.contour(grid1d, grid1d, density_values, cmap=\"plasma\",\n",
    "                # Because built-in heuristical levels cause animation noise:\n",
    "                levels=np.linspace(1e-4, 1/np.sqrt(la.det(2*np.pi*C)), 11))\n",
    "\n",
    "    # See exc. below\n",
    "    if seed and 'sample_GM' in globals():\n",
    "        plt.scatter(*sample_GM(mu, C=C, N=100, rng=seed))\n",
    "\n",
    "    plt.axis('equal');\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fdc4782",
   "metadata": {},
   "source": [
    "Note that the code defines the covariance `cv_xy` from the input ***correlation*** `corr`.\n",
    "This is a coefficient (number),\n",
    "defined for any two random variables $X$ and $Y$ (not necessarily Gaussian) as\n",
    "$$ \\rho[X,Y]=\\frac{\\mathbb{Cov}[X,Y]}{\\sigma_x \\sigma_y} \\,.$$\n",
    "Equivalently, it is the covariance between the *standardized* variables,\n",
    "i.e. $\\rho[X,Y] = \\mathbb{Cov}[X / \\sigma_x, Y / \\sigma_y]$.\n",
    "It quantifies (defines) the ***linear dependence*** between $X$ and $Y$,\n",
    "as illustrated by the following exercises.\n",
    "\n",
    "**Exc – Correlation influence:** How do the contours look? Try to understand why. Cases:\n",
    "\n",
    "- (a) correlation=0.\n",
    "- (b) correlation=0.99.\n",
    "- (c) correlation=0.5. (Note that we've used `plt.axis('equal')`).\n",
    "- (d) correlation=0.5, but with non-equal variances.\n",
    "\n",
    "Finally (optional): why does the code \"crash\" when `corr = +/- 1`? Is this a good or a bad thing?  \n",
    "\n",
    "More generally, it can be shown that $\\rho^2$ is the proportion of the variance of $Y$\n",
    "captured/explained by a simple linear regression from $X$.\n",
    "\n",
    "<a name=\"Exc-–-correlation-extremes\"></a>\n",
    "\n",
    "#### Exc (optional) – Correlation extremes\n",
    "\n",
    "Show that\n",
    "\n",
    "- (a) $\\rho[X,Y] = 0$ if $X$ and $Y$ are independent.\n",
    "- (b) $\\rho = 1$ if $Y = a X$ for some $a > 0$.\n",
    "- (c) $\\rho = -1$ if $Y = a X$ for some $a < 0$.\n",
    "\n",
    "Otherwise, it can be shown by Cauchy-Swartz, that $-1\\leq \\rho \\leq 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fc58ae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('Correlation extremes', 'a')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d54c6ccd",
   "metadata": {},
   "source": [
    "**Exc Correlation game:** [Play](http://guessthecorrelation.com/) until you get a score (gold coins) of 5 or more.  \n",
    "\n",
    "**Exc – Correlation disambiguation:**\n",
    "\n",
    "- What's the difference between correlation and covariance (in a single sentence)?\n",
    "- What's the difference between non-zero (C) correlation (or covariance) and (D) dependence?\n",
    "  *Hint: consider this [image](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#/media/File:Correlation_examples2.svg).*  \n",
    "  - Does $C \\Rightarrow D$ or the converse?  \n",
    "  - What about the negation, $\\neg D \\Rightarrow \\neg C$, or its converse?*  \n",
    "  - What about the (jointly) Gaussian case?\n",
    "- Does correlation (or dependence) imply causation?\n",
    "- Suppose $x$ and $y$ have non-zero correlation, but neither one causes the other.\n",
    "  Does information about $y$ give you information about $x$?\n",
    "\n",
    "<a name=\"Exc-–-linear-algebra-with-random-variables\"></a>\n",
    "\n",
    "#### Exc – linear algebra with random variables\n",
    "\n",
    "- (a) Prove the linearity of the expectation operator:\n",
    "  $\\Expect[a X + Y] = a \\Expect[X] + \\Expect[Y]$.\n",
    "- (b) Thereby, show that $\\mathbb{Var}[ a  X + Y ] = a^2 \\mathbb{Var} [X] + \\mathbb{Var} [Y]$\n",
    "  if $X$ and $Y$ are independent.  \n",
    "- (c) Similarly, prove:\n",
    "  $\\mathbb{Cov}[ \\vect{A} \\, \\vect{X} + \\vect{Y} ] = \\mat{A} \\, \\mathbb{Cov} [\\vect{X}] \\, \\mat{A}\\tr + \\mathbb{Cov}[\\vect{Y}]$ if $\\vect{X}$ and $\\vect{Y}$ are independent.\n",
    "- (d – *optional*) If $X$ and $Y$ are Gaussian, then so is $X + Y$.\n",
    "  Proof in the [next tutorial](T3%20-%20Bayesian%20inference.ipynb#Exc-–-BR-LG1). Meanwhile watch the [`3blue1brown` video](https://www.youtube.com/watch?v=d_qvLDhkg00&t=266s&ab_channel=3Blue1Brown).\n",
    "- (e) Let $\\vect{Z} \\sim \\NormDist(\\vect{0}, \\I)$,  where $\\I$ is the identity matrix.\n",
    "  Show that each component, $Z_i$, is independent of all others."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa3180ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('RV linear algebra', 'a')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97fcb678",
   "metadata": {},
   "source": [
    "#### Exc – Gaussian (multivariate) sampling\n",
    "\n",
    "Pseudo-random samples of $Z_i\\sim \\NormDist(0, 1)$\n",
    "can be generated on any modern computer using one of [these algorithms](https://en.wikipedia.org/wiki/Normal_distribution#Computational_methods).\n",
    "As shown above, the \"transformation\" $\\X = \\mat{L} \\vect{Z} + \\mu$\n",
    "yields $\\vect{X} \\sim \\NormDist(\\mu, \\mat{C})$,\n",
    "which can be used to sample with a desired mean and covariance, $\\mat{C} = \\mat{L}^{}\\mat{L}^T$.\n",
    "Indeed, the code below samples $N$ realizations of $\\X$,\n",
    "and assembles them as columns in an \"ensemble matrix\", $\\E$.\n",
    "But `for` loops are slow in Python (and Matlab).\n",
    "Replace it with something akin to `E = mu + L@Z`.\n",
    "*Hint: this snippet will fail because it's trying to add a vector to a matrix.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ee391dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_GM(mu=0, L=None, C=None, N=1, reg=0, rng=rnd):\n",
    "    # Seed random number generator\n",
    "    if isinstance(rng, int):\n",
    "        rng = rnd.default_rng(seed=rng)\n",
    "\n",
    "    # Compute L from C (if needed)\n",
    "    if L is None:\n",
    "        from numpy.linalg import cholesky\n",
    "        if reg:\n",
    "            C = C + reg * np.eye(len(C))\n",
    "        L = cholesky(C)\n",
    "\n",
    "    d = len(L) # len (number of dims) of x\n",
    "    Z = rng.standard_normal((N, d)).T\n",
    "\n",
    "    # Ensure mu is 1d\n",
    "    if np.isscalar(mu):\n",
    "        mu = mu * np.ones(d)\n",
    "\n",
    "    # Using a loop (\"slow\"):\n",
    "    E = np.zeros((d, N))\n",
    "    for n in range(N):\n",
    "        E[:, n] = mu + L @ Z[:, n]\n",
    "\n",
    "    return E"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8901db5a",
   "metadata": {},
   "source": [
    "Go back up to the interactive illustration of the 2D Gaussian distribution re-run its cell to check (eyeball measure) your implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bb7b588",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('Broadcasting')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f468c4fa",
   "metadata": {},
   "source": [
    "**Exc (optional) – Gaussian ubiquity:** Why are we so fond of the Gaussian assumption?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e43ac120",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show_answer('Why Gaussian')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ace12e0d",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "The Normal/Gaussian distribution is bell-shaped.\n",
    "Its parameters are the mean and the variance.\n",
    "In the multivariate case, the mean is a vector,\n",
    "while the second parameter becomes a covariance *matrix*,\n",
    "whose off-diagonal elements represent scaled correlation factors,\n",
    "which measure *linear* dependence.\n",
    "\n",
    "### Next: [T3 - Bayesian inference](T3%20-%20Bayesian%20inference.ipynb)\n",
    "\n",
    "<a name=\"References\"></a>\n",
    "\n",
    "### References\n",
    "\n",
    "<!--\n",
    "@book{laplace1820theorie,\n",
    "title={Th{\\'e}orie analytique des probabilit{\\'e}s},\n",
    "author={de Laplace, Pierre Simon},\n",
    "volume={7},\n",
    "year={1820},\n",
    "publisher={Courcier}\n",
    "}\n",
    "\n",
    "@book{gauss1877theoria,\n",
    "title={Theoria motus corporum coelestium in sectionibus conicis solem ambientium},\n",
    "author={Gauss, Carl Friedrich},\n",
    "volume={7},\n",
    "year={1877},\n",
    "publisher={FA Perthes}\n",
    "}\n",
    "-->\n",
    "\n",
    "- **Laplace (1812)**: P. S. Laplace, \"Théorie Analytique des Probabilités\", 1812.\n",
    "- **Gauss (1809)**: Gauss, C. F. (1809). *Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium*. Specifically, Book II, Section 3, Art. 177-179, where he presents the method of least squares (which will be very relevant to us) and its probabilistic justification based on the normal distribution of errors."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "jupytext": {
   "formats": "ipynb,nb_mirrors//py:light,nb_mirrors//md"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
